The inverse Laplace transform as the ultimate 
tool for transverse mass spectra 

Ekkard Schnedermann 
Physics Department, Brookhaven National Laboratory 
Upton, New York 11973, USA 

^ 20 Jan 1994 

P-h Abstract 

CN New high statistics data from the second generation of ultrarelativistic 

heavy-ion experiments open up new possibilities in terms of data analy- 
sis. To fully utilize the potential we propose to analyze the mj_-spectra 
of hadrons using the inverse Laplace transform. The problems with 

CN its inherent ill-definedness can be overcome and several applications 

in other fields like biology, chemistry or optics have already shown its 
feasability. Moreover the method also promises to deliver upper bounds 
on the total information content of the spectra, which is of big impor- 

+ j > tance for all other means of analysis. Here we compute several Laplace 

f q inversions from different thermal scenarios both analytically and numer- 

^ ically to test the efficiency of the method. Especially the case of a two 

C component structure, related to a possible first order phase transition 

to a quark gluon plasma, is closer investigated and it is shown that at 
least a signal to noise ratio of 10 4 is necessary to resolve two individual 
components. 

1 Introduction 

The hadronic transverse momentum spectra are one of the easiest accessible 
observables in relativistic heavy-ion collisions and are currently the main source 
of information about the dynamical evolution of the collision zone. With the 
arrival of new high statistics data from the second generation of heavy-ion 
experiments at CERN [2] and the planned high statistics experiments at RHIC 
(Relativistic Heavy Ion Collider) [3] new prospects for their analysis open up. 
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Customarily the analytical approach consists of performing a least squares 
fit of an exponential function to the measured mj_-spectra in order to deter- 
mine the global slope, which might be directly or indirectly connected to the 
temperature of an equilibrated hadronic system. It has already shown its lim- 
itations in the sense that the current data cannot be described by a single 
exponential over the whole range in m±_ [4, 5]. Fitting instead a sum of two or 
even more exponentials (along with their normalizations) to the data becomes 
an increasingly ill-defined procedure as the possible errors become strongly 
correlated. 

A straightforward extension of the method is to analyze the local slopes as 
a function of m±_ instead of the global slope. It can be performed for better 
theoretical insight into the spectral shapes from resonance decays or transverse 
flows, e.g. [6]. However, it seems to be impractical experimentally, because the 
approximation of the exponential slopes by finite differences between the mea- 
sured data points will lead to huge error bars. Moreover, at least in our opinion 
as opposed to [7], the local slopes plotted over m±_ have no obvious physical 
meaning, thus making it difficult to interpret the slopes without resorting to 
an actual model calculation. 

We want to propose the inverse Laplace transform for the analysis of almost 
exponentially decaying spectra like the mj_-spectra. It performs the analysis 
of a function in terms of decaying exponentials quite analogous to the Fourier 
transform, which employs oscillating exponentials. Though its use has been 
suggested several times in the past, the Laplace transform has not yet matched 
the tremendous success of the Fourier transform in data analysis; the numerous 
literature, e.g. [9], is rather focused on solving differential equations. Partly 
that is because the measurement has to be precise over several orders of mag- 
nitude, in analogy to several oscillations needed for the Fourier transform. But 
the main reason is that the construction of the inverse Laplace transform for 
real data is an inherently ill-defined or ill-conditioned problem [9, 11]. This 
problem has been solved by means of an eigenfunction or singular function 
analysis [12, 13] and many applications from diverse as medicine, 

biology, chemistry and optics have been investigated [14], amonst them the re- 
trieval of the temperature distribution inside the human body using microwave 
radiation. 

After showing the problems with naive fit procedures we want to investigate 
in this paper the applicability of the inverse Laplace transform to ultrarela- 
tivistic heavy ion collisions. Assuming that a locally equilibrated system is 
produced in these collisions, particle spectra can be computed for a number of 
scenarios. Most of them can be inverted analytically to provide us with a set of 
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patterns to look for in the inversion of real data. The actual inversion of data 
by means of a singular function analysis is then performed for several test- 
cases, focussing especially on the requirements for resolving a two component 
structure in the spectra. 

2 Naive Methods 

As a practical exercise we perform fits by sums of exponentials with c 4 - and t{ 
as the free parameters 

n 

g^tim^) = c i exp(-m_L t,) (1) 

8 = 1 

to a test spectrum, which we construct from two thermal distributions with 
different temperatures T\ and T 2 integrated over rapidity (or fluid rapidity in 
a boost invariant situation): 

dn m_i_ frn±\ m±_ fm± 

gt es t{m L ) = oc -=-Ki [-=- J + -zr&i ~7f- 

m^dm^ ii V ii / ±2 V ±2 



where mj_ = Jmg + p]_ is the transverse mass of the particle. 
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Table 1: Fitting the spectrum (2) by a sum of exponentials (1) gives a 
stable fit in the case of two slopes with a good agreement (x 2 ~ 2 X 10~ 5 ) 
to the test function, which however does not reproduce the built in 
components at / = 0.9/200 MeV and / = 1.1/200 MeV. A three slope fit 
does not improve this situation but is already difficult to perform due 
to the many local minima in the six dimensional parameter space from 
which we are showing the three fits with % 2 ^ 10~ 6 



For a specific example we have chosen m = 140 MeV, T\ = 200/0.9 MeV 
and T 2 = 200/1.1 MeV so that the two components should appear at t\ = 
0.9/200 MeV" 1 and t 2 = 1.1/200 MeV" 1 . Note that we adjusted the relative 
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normalizations so that their Laplace components (12) are roughly equal in size, 
although the number of particles from the higher temperature dominate the 
ones from lower temperature by a factor of 1.56. The quality of the fits is then 
judged by summing the squared difference over 32 data points pi spaced at 
50 MeV intervals: 
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Table 2: Fits to the spectrum (2) by sums of exponentials (1) with the 
slopes t{ fixed and the normalizations c 4 - as the free parameters. The 
obvious instabilities in the solution grow with the number of parameters 
and render this method unusable for higher precision data. 



We performed two different fit procedures. In the first case, shown in 
Tab. 1, we performed a fit to the spectrum (2) with both the c 4 - and t{ as free 
parameters. Because of the many local minima involved this becomes a rather 
tricky fit for more than three exponentials so that in the second case, shown 
in Tab. 2, we fixed the slopes t{ at chosen values and only performed fits with 
the normalizations free parameters. 

From the study of the tables it becomes apparent that the fits are becom- 
ing more and more unstable when the number of parameters is increased and 
are thus meaningless. In no case were the fits able to capture the built-in two 
component nature of the spectrum, which we ascribe to the high sensitivity of 
this prescription to small "fluctuations" in combination with the fact that the 
Bessel functions are only approximately exponential. In the absence of noise 
and for pure exponentials there is apparently no problem, although the % 2 con- 
tours are elongated ellipsoids, thus showing very big correlated errors between 
the two slope parameters. In the following sections we will demonstrate that 
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the inverse Laplace transform by the method of singular functions is much 
better suited to the occurrence of deviations from exponential behaviour and 
noise in the data. 

3 Analytical Inversions 

The Laplace transform of the function fit) is given by 

oo 

g(p) = C{f} = j f(t) e"* dt (4) 
o 

and is a linear transformation. This has the effect that different Laplace com- 
ponents fit) = fi(t) with their corresponding spectra gi(p) will give a spec- 
trum which again is the sum of its components g(p) = 9i(p)- This seemingly 
trivial property is not valid for the analysis in terms of local slopes, so that the 
connection between a certain slope value and the underlying physics is rather 
hidden. 

In the case of an extended static source in local thermal equilibrium we 
would interpret the Laplace components fit) as the temperature distribution 
(actually I/T-distribution) of the source. For hadronic spectra we would ex- 
pect that they are concentrated around the freeze-out temperature, where the 
equilibrium between particles can no longer be kept up and they are streaming 
freely into the detectors. The appearance of a two component structure in 
the inverse fit) would hint towards the presence of a phase transition into the 
quark gluon plasma. The second peak, located at t = 1/T C , could originate 
from the hadrons which escape from the mixed phase during its long life time 
spent at the critical temperature T c , and travel to the detectors without further 
interactions in the surrounding hadron gas. 

To obtain the components fit) from the measured spectrum g(p) } or g(m±) } 
we have to perform the inverse Laplace transformation. It is given by Bromwich's 
integral in the complex plane 

c-H'oo 

f(t) = C- 1 {g} = — I g(p)e-*dp (5) 

Z7TZ J 

c—ico 

with c any real number so that g(p) is analytical for Rep > c . We will use this 
method, or a numerical implementation of it via the NAG-Library, to compute 
the inverse Laplace transforms for several spectral shapes which follow from 
theoretical scenarios. However, for measured data a different procedure has to 
be developed in section 4, because they are not available in the complex plane. 
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Nevertheless, inversions using eq. (5) can be performed for many theoret- 
ical scenarios, if the theoretically computed spectra have unique extensions 
into the complex plane. We give inverse transforms here for several kinds of 
(hadron) mj_-spectra, because we think they have the best chances of appli- 
cation, though also pj_-spectra and (dilepton) invariant mass spectra might be 
conceivable. We list several formula (with the help of [10]), which relate to 
transforming the spectrum shifted by the restmass (7), multiplication by m±_ 
(8), Boltzmann distributions (10), integrated thermal distributions (12) and 
Bose/Fermi-distributions (13). 

C-'igim^)} = /(f) (6) 

£ _1 {#(ra_L)} = exp(-tm )C~ 1 {g(m 1 _ - m )} (7) 
d 

C 1 {m L g{m i _)} = — C 1 {g(m L )} - ljm_ C 1 {g(m L )} (8) 

£- 1 {exp(-m 1 /T)} = 6(t - ±) (9) 

£- 1 {m 1 exp(-m 1 /r)} = S'(t - i) (10) 

f for < t < ^ 

C-^MmjT)} = U fort> i (11) 

{ V(^) 2 -i T 

f for < t < i 

C-^m^MmjT)} = \s(t-±)-^^ fort>i ( 12 ) 

c ~ 1{ — r7nT7 } = h^y-^'it-^) (13) 

exp(m_L/r) =p 1 £r[ v J 1 / 

The spectral shape of (12) involving the Bessel function K\ results from in- 
tegrating the thermal spectrum of a point source over rapidity to compute 
dn/dm 2 ^, or equivalently, to integrate over the fluid rapidity rj in a boost in- 
variant scenario to compute dn/dydm 2 ^ (e.g. 



4 Results 

In the case of real data, which are naturally given only on the real axis, we have 
to proceed differently, because inverting the Laplace transform for functions 
which are defined only on the real axis instead of at least a complex half plane is 
an ill-conditioned problem [11]. For data points with associated measurement 
errors the problem even gets worse. Historically there have been many naive 



6 



attempts to compute the inverse Laplace transform, e.g. by fitting the data 
with superpositions ol exponentials (e.g. [8]) 



n 



f(t) = Y, c ^(t-t t ), 



(14) 



or combinations ol exponentials and polynomials [15]. 

This is not a viable method because the inverse Laplace transform is not 
continuous if no further restrictions on the possible space of ordinary functions 
are imposed. This can be illustrated by the following Laplace transform [9]: 



which can be made arbitrarily small with u — > oo. One could thus approach 
with Cf the given data while / is not at all converging to a certain solution. 
This can be cured by expanding / in terms of a series of singular functions [11, 
12, 13]. This effectively restricts the available functional space for / by cutting 
out high frequency oscillations (=noise), which would otherwise drastically 
influence the solution. 

We follow the method given in [13], which assumes that / is roughly dis- 
tributed like some profile function P(t), which is in this case chosen to be the 
Gamma distribution with width f3: 



This profile function can be tailored to the specific problem at hand, even step 
functions can be chosen [12]. The width f3 can be taken, as we did, from pure 
prejudice about what the components should be. The choice can actually be 
checked, because if the window profile is chosen too small the inversion will try 
to fold the outside components into the window which will in turn lead to big 
oscillations and will depend strongly on the actual number of singular values 
used. However, ideally the width would have to be estimated from the data, 
e.g. by the cumulants method [15], which is however a rather complicated 
procedure in itself. 

In terms of data sampling we are not as diligently adjusting the spacing of 
the data points as it is done in [13] in order to retrieve maximal information 
content. We simply assume that the data consists of 32 points taken at 50 MeV 
intervals, though apparently exponential sampling [16] is more efficient, espe- 
cially when it is matched to the width of the profile function. However, our 



C{sinujt} 




(15) 




(16) 
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choices are not very far from optimal because our necessary signal to noise 
ratio seems to be in line with [13]. 

Because there is only a finite number of sampling points, the Laplace inver- 
sion by means of singular functions turns into a singular value decomposition 
for ordinary matrices, which can easily be performed numerically and which 
yields a descending series of singular values a,. By truncating that series at a 
certain a n , the part of the solution which leads to wild oscillations is filtered 
out and one obtains the approximate solution as a series of singular func- 
tions. The precision of the regularized solution, i.e. the noise to signal ratio, 
is reflected in the condition number a n ja\. Naturally, also more elaborate 
regularization schemes can be used [17]. 
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Figure 1: Numerical inversion of 
the spectral function as performed 
by the NAG-Routines. We assumed 
a spread in the thermal components 
of the size of 5%. The particle's 
mass has been fixed as the pion 
mass. 
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Figure 2: Reconstruction of two 
thermal pion distributions with 
T f = 180 MeV and T c = 220 MeV 
using 32 data points spaced at 
50 MeV and a very narrow profile 
function with 3 = 20. A clear 
seperation of the two components is 
seen for the case of 5 singular func- 
tions, corresponding to a noise level 
of a 5 /a t < 1.5 X 10~ 4 . 



Some general examples of inversions are shown in [13]. Here we will con- 
centrate on several cases which are more closely related to possible heavy-ion 
situation scenarios and we return to our initial motivation, which is to search 
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for the phase transition by looking for a two component structure. In Figs. 1 
and 2 we investigate the resolution limits of the inversion method for two ther- 
mal distributions which are very close together (2), i.e. Tj = 182 MeV for the 
freeze-out component and T c = 222 MeV for the critical temperature of the 
phase transition. Especially with respect to T c those values are meant only as 
moderately realistic examples which are modeled after the currently measured 
spectra. Note that the apparent temperatures do not have to be the actual 
local temperatures in the system since radial flows of average (v r ) 0.3 c give 
rise to a blueshift of about 30% [6]. 

First, in Fig. 1, we have used the N AG-Library to perform a Laplace inver- 
sion of a two-component thermal spectrum of pions, which has been extended 
to the complex plane, to have a benchmark for the following inversion of data 
on the real axis. In Fig. 2 we then perform a Laplace inversion from a singular 
value decomposition, chosing a very narrow profile function (ft = 20) centered 
at t = 1/200 MeV -1 . We see that we indeed can seperate both components 
much better than with the naive attempts section 2, though the reproduction 
of the shape of the actual inverse Laplace transform is rather marginal with 
this noise level. Nevertheless, this scenario shows the basic applicability of the 
method to actual heavy-ion data and a promising new way to analyze high 
statistics data. 

However, we note that the actual heavy-ion collision might produce apart 
from the purely thermal radiation also resonances which afterwards decay into 
pions [18]. Again to enable comparison with Fig. 4 we perform in Fig. 3 the 
inversion of an integrated thermal distribution (12) of pions at T = 200 MeV 
using a profile function centered at T with a width ft = 2, which is much 
broader than the actual distribution. We note that the reconstruction, while 
much broader than the actual inversion, is still narrower than the profile func- 
tion. Moreover we do not attribute any physical meaning to the fluctuations 
at low t, which might simply reflect an unfortunate choice of profile function 

P(t)- 

In Fig. 4 the additional pions from decays of a thermal population of p 
and to resonances (modeled after [18]) clearly gives rise to a second component 
which is well seperated from the first one. Its appearance at larger values of 
t is linked to the much lower apparent temperatures T app = 1/t of resonance 
decay spectra when compared to the actual temperature of the system [18, 6]. 
The inversion requires at least five singular functions (shown is seven) and a 
noise to signal ratio of at least a 5 /cti £ 5 X 10~ 3 (ctr/cti ^2 X 10~ 4 ). 

Apparently other processes, e.g. resonance decays, can also simulate a sec- 
ond temperature component, which is a drawback to the proposed method. 
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Figure 3: Laplace inversion of a 
thermal pion m_|_-spectrum using 32 
datapoints spaced at 50 MeV and 
using the first seven singular values 
a, = 0.71, 0.20, 0.059, 0.015, 0.0036, 
0.00081 and 0.00017. The profile 
P(t) was chosen very wide {(i = 2). 
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Figure 4: The same Laplace inver- 
sion as in Fig. 3 now for a pion spec- 
trum additionally including p and 
uj decays shows a second compo- 
nent at higher t. This reconstruc- 
tion corresponds to a noise level of 
ar/«i £ 2 x 10" 4 . 



Nevertheless, there are two possibilities around this obstacle, because the de- 
cay contributions necessarily depend on the kind of particle measured and the 
produced resonances. By selecting particle types which are less sensitive to the 
resonance contribution than the pions, e.g. kaons, and comparing them among 
one another, one could eliminate large parts of that background, because only 
the thermal components should look similar across all particle species. On 
the other hand, if there are very many resonance channels, the multitude of 
channels might rather form a continuum of Laplace components, above which 
certain peaking thermal components at Tj and T c might still be identifiable. 

5 Summary and Conclusion 

We have investigated the applicability of the inverse Laplace transform on 
heavy-ion data and have shown that naive methods like fitting multiple ex- 
ponentials are incapable to analyze high statistics data. Using the method 
of singular function decomposition [13] we have shown how a two component 
structure in the produced hadron spectra can be identified and that the re- 
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quirements for the data in terms of noise level are within reach of current or 
planned experiments. We also would like to argue that the issue of systemat- 
ical errors becomes less severe because the method projects some errors into 
certain regions so that they are eliminated or identifyable. However, only a 
simulation of the actual experiments could assure that. 

However we have also shown how resonance decays can cloud the picture by 
generating a second component at about twice the thermal slope. Performing 
an analysis with kaons or baryons should reduce this component and also 
change its shape so that the whole body of data should give clear answer. 

We have tried a Laplace inversion with tt + and tt~ data from the ISR 
[19] which have statistical errors of about 1-2%. However, we could not find 
any indication of a two component structure, which might either indicate the 
absence of a (first order) phase transition or insufficient data to decide the 
question. 

The hadron spectra are the most promising for an analysis because of their 
high statistics. The inversion of the lepton spectra might be even more inter- 
esting, e.g. one would look for a single component at T c above a rather smooth 
background, but the necessary precision seems to be far out of reach. 

Further improvements of the singular function method, e.g. exponential 
sampling, choice of profile function, can conceivably improve upon the effi- 
ciency of the method, but they will require a close relation to the performed 
experiment. 
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